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We study the extended Bose-Hubbard model describing an ultracold gas of dipolar molecules 
in an optical lattice, taking into account all on-site and nearest-neighbor interactions, including 
occupation-dependent tunneling and pair tunneling terms. Using exact diagonalization and the 
multiscale entanglement renormalization ansatz, we show that these terms can destroy insulating 
phases and lead to novel quantum phases. These considerable changes of the phase diagram have 
to be taken into account in upcoming experiments with dipolar molecules. 
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Trapping and manipulating ultracold gases in optical 
lattices has allowed the realization of many-body physics 
in a controlled environment. For atoms interacting via 
contact interaction, a quantum phase transition from a 
superfluid (SF) to a Mott insulator (MI) has been pre- 
dicted and observed In the simplest case, these sys- 
tems can be theoretically described by the Bose-Hubbard 
(BH) model, which has two parameters: a tunneling J 
and an on-site interaction ?7 @, 0| • A natural extension 
of the Bose-Hubbard model comes from including long- 
range interactions between particles. Experiments on ul- 
tracold polar molecules have renewed interest in extended 
Bose-Hubbard models which can model such systems in 
optical lattices • Because of the strong electric dipole 
moment of polar molecules, long-range interactions play 
a crucial role in the collective behavior of the system, 
leading to the appearance of states with long-range or- 
der, like various structured insulating states, supersolids, 
Wigner crystals, pair-supersolids, etc. 



In this Letter, we study the ground-state of dipolar 
molecules in a 2D square optical lattice with a harmonic 
trapping along the polarization direction of the dipoles. 
We derive a modified BH model which includes addi- 
tional occupation-dependent nearest-neighbor (NN) hop- 
ping processes arising from long-range dipolar interac- 
tions in the lowest Bloch band. Usually, interaction- 
induced hopping terms are neglected when discussing 
dipolar bosonic molecules. In this Letter, we show that 
these terms considerably change the physics of dipolar 
soft-core bosons. Soft-core bosons in square and one- 
dimensional lattices have been discussed in the litera- 
ture within the extended Hubbard model, focusing on 
the presence of stable supersolidity 
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In the 

usual case with only NN interaction, at sufficient dipo- 
lar strength, the ground states at half- and unit-filling 
are checkerboard (CB) insulating states. Using exact di- 
agonalization (ED) and multiscale entanglement renor- 
malization ansatz (MERA) , we solve the one-dimensional 
extended Hubbard model including the novel occupation- 
dependent NN hopping processes. We find that with in- 



creasing dipolar interaction, the system enters from the 
CB phases to a novel state which has a one-particle su- 
perfluid (SF) and pair-superfluid (PSF) properties. Par- 
ticularly we find a region where both of them coexists 
with the SF order parameter has alternating sign at con- 
secutive sites. 

Our system consists of dipolar bosons polarized by an 
external electric field along the z direction and confined 
in a square optical lattice. The corresponding Hamilto- 

nian reads H = /d^r ^'t(r) -^V^ -I- Matt(^) "^{r) + 

i//d3rd3r'*t(r)^ft(r')V(r - r')*(r)*(r'), where 
^'t(r) (^'(r)) are the bosonic creation (annihilation) field 



operators. Vutt{r) = Vb [sin 
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sin^ ^y\+mnlz^ /2 



is an external lattice potential of lattice depth Vq, gen- 
erated by a laser field of wave-length A, with 17^ charac- 
terizing the external harmonic potential in z direction. 
The dipole-dipole interaction is denoted by V(r). By ex- 
panding the field operator "^{r) = ^ - Wi{x, y)e~'^^ 1"^ en 
in lowest Bloch-band Wannier- functions yVi(a;, y), and by 
restricting ourselves to on-site and NN terms, we arrive 
at the extended BH model 



= - J ^ aja, + ^ n,(ri, - 1) -t- ^ 



niUj 



(1) 



(a|) annihilates (creates) a particle on lattice 



- is the corresponding density operator, J 



where a. 
site i, fii = a^ai 
the standard tunneling coefficient, U the on-site interac- 
tion, and V the NN interaction, arising from a truncation 
of the dipolar interactions to the dominating term. Dipo- 
lar interactions lead to two novel terms in Eq. ([T]): The 
term proportional to T describes one-particle tunneling 
to a neighboring site induced by the occupation of that 
site, and the term propo rtional to P is responsible for 
NN pair tunneling [l9l-[21|. 

The matrix elements U , V, T, and P are given by a 
sum of dipolar and (5-like contact interactions, V(r— r') = 



FIG. 1. Dependence of the dipolar part (subscript D) of 
U, V, T, and P on the lattice flattening k for lattice depth 
Vo = 6Er and 7 = 52. 
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We measure all 



lengths in units of the laser wave length A and all en- 
ergies in recoil energies Ej^ = 2TT'^h'^ / (mX^), where m 
is the bosonic mass. Additionally, we define the lattice 
flattening k = hD,z/2E]^ as well as the dimensionless cou- 
pling constants describing contact and dipolar interac- 
tion, g = Idn'^as/X and 7 = mcP /{h^eaX) (where is 
the s-wave scattering length, Eq is the vacuum permittiv- 
ity, and d is the electric dipole moment of the bosons). 

For concreteness, we consider an ultracold gas of dipo- 
lar molecules confined in a optical lattice with lattice 
depth Vb = QEji, mass m = 220a.m.u and A = 790 nm 
[221 . We also assume that the s-wave scattering length 
of the molecules, ~ lOOoo. For these parameters, 
g w 1.06 is approximately constant. We consider dipole 
moments d up to ^ 3D (7 up to ~ 470), which can be 
achievable for molecules like bosonic RbCs, KLi Q etc. 
To illustrate the relative strengths of different parame- 
ters, in Fig. [U we compare for 7 = 52 the tunneling J 
with the dipolar contribution (subscript D) to the pa- 
rameters t/, V, T and P. For the parameters chosen, Td 
and Pq are 1 orders of magnitude smaller than Vd where 
as Uu/Td can be tuned by changing k. On the other 
hand, Td can dominate over J for large 7. In addition, 
T and J can have opposite sign as seen in Fig. [1] For 
concreteness, we choose the lattice parameter k w 1.95, 
making (additionally to J) the on-site interaction U al- 
most independent of the dipole moment (Uu ~ 0). In this 
case, for large enough 7, we expect that with increasing d 
the parameters V, T and P determine the system prop- 
erties. For clarity, we restrict ourselves to a ID chain of 
N lattice sites with periodic boundary conditions. 

To get a first understanding of the system, we find 
the ground state \^Q{d)) as a function of d by exact di- 
agonalization (ED) of a half-filled system with = 8 
sites. We also present results for iV = 12 and 16 to check 
for dependence on system size. Without the occupation- 
dependent tunneling terms T and P, we observe the usual 
scenario with only two phases, a single-particle SF and 
a CB phase. The transition happens at d « 0.4 D. It is 
marked by an increase of the contribution of the checker- 
board states to the ground state to almost 100% [inset 
of Fig. El^a)]. Also, the one-particle correlation function 
cj)i = almost vanishes, indicating the transi- 
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FIG. 2. We plot various properties of the exact ground state of 
a half-filled system as a function of dipole moment d. Fig. (a) 
shows the contribution of the CB states to the ground state 
of the system for N = 8. In Fig. (b) we plot the one-particle 
and two- particle correlation functions (j>i and The dotted 
line shows (pi when we neglect the terms T and P. When 
T,P ^ 0, the solid line and dash-dotted line shows and 
respectively as a function of dipole moment d. In Fig. (c) we 
plot the fidelity susceptibility x{d) for the half-filled system 
for different system sizes. In Fig. (d) we have shown the 
structure-factor S at different ordering wave vectors for the 
half-filled system with 16 sites. 



tion to an insulating state. In the half-filled system, the 
transition occurs because for large enough V the parti- 
cles can decrease their energy by avoiding every second 
site. If we neglect T and P, the situation will not change 
by further increasing d (dotted lines in Fig. [5]), since 
this only increases V even more. However, the situa- 
tion changes significantly when we take into account the 
density-induced tunneling T and the pair tunneling P. In 
this case, for d w 1.1 D, a second phase transition occurs, 
destroying the CB order [solid lines in Fig. [2Ia)[. Pre- 
vious studies have completely neglected such a possible 
destruction of CB order at large d. At the transition, the 
contribution of the CB state to the ground state decreases 
rapidly, and the one-particle as well as the two-particle 
NN correlation function $i — [dashed- 
dotted line in Fig. [2l^b)[ attain finite positive values, in- 
dicating that the new phase shows single-particle as well 
as pair superfluidity. In this region we also find that the 



long-ranged correlation function (aja.;) for \i — j\ < 6 
decays slowly with alternation sign for consecutive sites. 
This suggests appearance of antiferromagnetic like order 
due to the positive hopping T resulting in the conden- 
sation of bosons at the edge of the Brillouin zone. We 
also looked into the relative effect of T and P on the PSF 
state. We found that PSF is generated due to the cor- 
related tunneling term T (in interplay with the nearest- 
neighbor interaction V). 

For even larger electric moments, a third phase transi- 
tion happens, where 1/)^ changes sign. Another signature 
of this transition is a rapid growth of 4>i. Since this 
quantity measures fluctuations of bosonic pairs, this is 
a signature of a novel pair-superfluid (PSF) phase. The 
appearance of pair superfluidity has previously been pre- 



dieted in bilayer dipolar systems where the particles are 
bound by an attractive interaction between the layers 
Though in bilayer systems, the state is a true 
molecular superfluid as $i ^ 0, whereas 0^ = iden- 
tically. In the present system, in spite of the particles 
interacting repulsively, the pairs are created due to the 
occupation-dependent tunneling terms in Eq. ([T]) (similar 

to [ii). 

To confirm that all these transitions are indeed quan- 
tum phase transitions, we calculated - for different 
chain lengths N - the ground-state fidelity suscepti- 

where T{d, 6) = 



T = P = 



S=0 



bihty mm X{d) = 

\{'4'Q{d)\'ipQ{d + <^))|. Peaks in x a-re efficient indicators of 
quantum phase transitions. In Fig. [JJc), we present xi^) 
for different chain sizes. There are three clear peaks at 
the quantum phase transitions found from the correlation 
functions [as presented in Fig.[2][b)]. The positions of the 
transition points (TPs) do not significantly depend on the 
number of sites, especially for the 1st and 3rd TP. The 
middle peak in Fig. [2Kc) refers to the transitions from 
checkerboard to antiferromagnetic superfluidity. More- 
over, the magnitude of the fidelity susceptibility at all 
TPs increases with chain length, which suggests that the 
transitions will survive in the thermodynamic limit. 

More insight into the properties of the observed phase 
comes from the static structure factor, which is defined 



with 



as S{q) = ^T.lk^i^''^'-''H{njnk) - {nj){nu)), 
q = 2'Kra/N , Q < m < N ~1 integer. A peak in the struc- 
ture factor at finite momentum points towards presence 
of periodic density modulation in the systems. In Fig. 
[2](d), we present S{q) for a half-filled system with = 16 
sites. In the CB phase (between the 1st and 2nd TP), 
the dominant peak of S{q) is at g = tt, and its magnitude 
is almost independent of system size. Above the 3rd TP, 
the system is in a phase where (pi has an inverted sign 
and $i is large. This means that states where bosons oc- 
cur in pairs dominate (their contribution to the ground 
state is about 95%). Since, due to the dipolar interac- 
tions, boson pairs do not occupy neighboring sites, the 
system has some local structure, leading to a predomi- 
nant structure- factor peak at q = ^. The intermediate 
phase (between the 2nd and 3rd TP) has interesting prop- 
erties: the ground state of the finite system deforms its 
structure stepwise, changing the dominant q from tt to 
7r/2 by one quantum Ag = 2tt/N at a time. For N = 16, 
this leads to three changes in the dominant q. Since in 
an infinite system q can take every value between and 
27r, we expect in large chains a continuous change from 
the CB with q = tt to the two-particle SF with q — 7r/2. 

Finally, we analyze the influence of the additional 
terms T and P on the grand-canonical phase diagram, 
where the particle number is not conserved. For this, we 
add a chemical potential term hi to Hamiltonian 

(U)). In Fig. [31 we present the phase diagram as well as 
the average number of particles per site for ED calcula- 
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FIG. 3. ED phase diagram without (left) and with (right) 
taking into account T and P. The color denotes the super- 
fluidity fractions, and Neglecting T and P, for large 
enough d and jj, the system is always in an insulating phase 
and the average number of particles is a multiple of 1/2. CB 
(CB2) denotes a checkerboard phase where sites with and 
1 (2) particles alternate. Including the new terms, the insu- 
lating phases vanish for large enough d, and a PSF appears. 
We truncate the Hilbert space at a maximal occupation num- 
ber of 4 particles per site. We exclude data points where the 
occupation number becomes too high (white region). 
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FIG. 4. Results for the BH Hamiltonian ([T| in a chain with 
A*' = 128 using MERA with m = 8, revealing checkerboard 
(CB and CB2) order, as well as superfluid (SF and SF') and 
pair-superfluid phases (PSF). (a) Mean occupation, (b) mean 
SF order parameter, (c) mean PSF order parameter, and (d) 
mean NN density correlations. 



tions of 4 sites with occupation truncated at 4 particles 
per site. When the additional terms T and P are large, 
they destroy the CB phase, making place for a PSF. 

To get a more detailed analysis of larger systems 
than tractable in ED, we have performed a Multi-Scale- 
Entanglement-Renormalization-Ansatz (MERA) (irl - Egj 
computation of the phase diagram The MERA is a quasi- 
exact variational method that consists in postulating 
a tensor-network structure for the low-energy states of 
Hamiltonian ([1]), which in particular yields especially 
good results in critical phases, where other methods such 
as DMRG are very costly [23,[28|. 

The results are presented in Figs. [DJa-d), where we 
show, averaged over the chain, the occupation (n^), the 
SF order parameter (a^), the PSF order parameter (0^0^), 
and NN density-^density correlations (71,71^+1). The phase 
diagram extracted from these observables is sketched in 
Fig. [D(a). At low d, there is a single-particle SF, which 
gives way to CB phases for d > fi. Increasing d, the sys- 
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tern undergoes a transition to a SF phase, where initially 
for a range of w 0.2D one-particle superfluidity domi- 
nates (similar to the ED results), and afterwards pair 
superfluidity. At low we find a phase (SF') which has 
additionally to SF order (i.e., a finite {ui)) small nearest- 
neighbor density-density correlations. Hence, it has a 
local structure where sites with high and low occupation 
alternate. We checked that this phase is not due to phase 
separation. The novel aspect of this is that in the usual 
extended BH model with soft-core interactions stable su- 
persolidity appears only at the particle-doped region of 
the CB phase (TtI . [18| . For higher fi and d ~ 1 , we get 
a CB of two particles(CB2 phase) in the filled site. This 
behavior is a result of having low U so that it is ener- 
getically favorable than having one particle at each site. 
As already indicated by ED, the new terms T and P 
destroy CB order in favor of PSF phases, meaning that 
these terms cannot be neglected. We also checked at few 
points in the phase space of the PSF region to look for 
the sign of the SF order parameter as a function of lat- 
tice sites and we found the alternating sign as seen in ED 
calculations. 

To make better contact with experiment, we examine 
the disappearance of the CB pattern when the long-range 
part of the full dipolar interactions is taken into account, 
i.e., we replace the NN term in Hamiltionian ([1]) with 
ji^f^if^j- Using ED at half-filling for N = 16, 
we find that qualitatively the phase diagram does not 
change much with respect to our previous calculations 
with the simplified Hamiltonian ([T]) [compare Fig. [5]Jb)]: 
When the occupation-induced tunneling terms T and P 
are neglected, the CB phase remains stable for arbitrarily 
large d [Fig. [51(a)]. In contrast, when taking into account 
the tunneling terms T and P it disappears, making way 
for a PSF phase [Fig.[5{b)]. This happens even at smaller 
d than when truncating the interactions at NNs. Namely, 
the PSF phase appears for d > 0.7D. We also note that 
in Fig. Eljb), there is a kink in around d ^ 0.5D. 
This kink corresponds to the appearance of a crystal 
like phase with modulation |....200100200100.... >. A de- 
tailed discussion of this phase is outside the scope of this 
paper. We have further checked that counter-intuitively 
PSF arises predominantly due to correlated tunneling T. 
Without this term PSF phase can not be reached for rea- 
sonable electric moments. We also note that for very low 
dipolar strength $i has a small nonzero value. As seen 
in Figs.[2]^a), (b) and[5fa), (b), a small but finite $i is 
present as c? — > irrespective of the presence of T and P. 
This can be traced back to second-order processes due to 
J which can also give rise to pair correlations with small 
magnitude. 

In summary, we showed - based on ED and MERA 
- that commonly neglected terms in the extended BH 
model for dipolar molecules in optical lattices can 
lead to interesting new phenomena. We showed for 
a particular choice of optical-lattice parameters that 
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FIG. 5. The one-particle and two-particle correlation func- 
tions (j>i (solid line) and $i (dotted line) as a functions of 
dipole moment d when the full dipolar interactions are taken 
into account corresponds qualitatively to the calculations 
truncated at NNs (Fig. [2|). The large and negative (pi when 
terms T and P are taken into account indicate the break down 
of the CB phase to a PSF. Calculations for ED at half-filling 
with N = 16. 

occupation-dependent tunneling and pair tunneling (in- 
duced by long-range dipolar interactions) destroy insu- 
lating checkerboard phases for large enough electric mo- 
ments d, leading to a novel pair-SF phase. MERA re- 
sults suggest also that a supersolid phase could appear 
for 1/2 filling even in the hole-doped case. Any presence 
of additional weak trapping potential can result in shell 
like structures seen in usual BH model as long as local- 
density approximation is valid. We note that, as our nu- 
merical calculations in carried out in one dimension, the 
various superfluid correlations decay in a power law with 
distance. In this sense, the superfluid phases mention 
here will show quasi-long-range order in infinite systems. 
Our calculations are done for parameters experimentally 
achievable in the near future, and the changes to the 
phase diagram have to be taken into account in the inter- 
pretation of future experiments with dipolar molecules. 

This Letter was supported by the EU STREP NAME- 
QUAM, IP AQUTE, ERC Grant QUAGATUA, Spanish 
MICINN (FIS2008-00784 and Consohder QOIT), Caixa 
Manresa, and Marie Curie project FP7-PEOPLE-2010- 
IIF "ENGAGES" 273524, and AAII-Hubbard. T.S. ac- 
knowledges hospitality from ICFO. 



[1] M. Greiner et. al.. Nature 415, 39-44 (2002). 

[2] M.P.A. Fisher, P.B. Weichman, G. Grinstein, D.S. 

Fisher, Phys. Rev. B 40, 546 (1989). 
[3] D. Jaksch et. al., Phys. Rev. Lett. 81, 3108 (1998). 
[4] K.-K. Ni et. al.. Science 322, 231 (2008). 
[5] S. Ospelkaus et. al., Faraday Discuss. 142, 351 (2009). 
[6] K. Aikawa et. al.. New J. Phys. 11, 055035 (2009) . 
[7] J. Deiglmayr, Phys. Rev. Lett. 101, 133004 (2008). 
[8] A.C. Voigt et. al., Phys. Rev. Lett. 102, 020405 (2009). 
[9] T. Lahaye, C. Menotti, L. Santos, M. Lewenstein, T. 

Pfau, Rep. Prog. Phys. 72, 126401 (2009). 
[10] K. Goral, L. Santos, and M. Lewenstein, Phys. Rev. Lett. 

88, 170406 (2002). 
[11] B. Capogrosso-Sansone et. al., Phys. Rev. Lett. 104, 



5 



125301 (2010). 

[12] G. Pupillo et. al., Phys. Rev. Lett. 100, 050402, (2008). 
[13] D.W. Wang, M.D. Lukin, and E. Dernier, Phys. Rev. 

Lett. 97,180413 (2006). 
[14[ A. Arguelles, L. Santos, Phys. Rev. A 75, 053613(2007). 
[15| C. Trefzger, C. Menotti, M. Lewenstein, Phys. Rev. Lett. 

103, 035304 (2009). 
[16[ M. Eckholt and J. J. Garcia-Ripoll, Phys. Rev. A 77, 

063603 (2008). 

[17| P. Sengupta et. al., Phys. Rev. Lett. 94, 207202 (2005). 
[18| G.G. Batrouni, F. Hebert, R.T. Scalettar, Phys. Rev. 

Lett. 97, 087209 (2006). 
[19[ A. V. Gorshkov et. al. Nature Physics 6, 289 (2010). 
[20[ G. Mazzarella, S. M. Giampaolo, and F. lUuminati, Phys. 

Rev. A 73, 013625 (2006). 
[21[ P. Bader and U. R. Fischer, Phys. Rev. Lett. 103, 060402 



(2009). 

[22[ S. Kotochigova, E. Tiesinga, Phys. Rev. A 73, 041405(R) 

(2006) . 

[23[ O. Dutta et. al.. New J. Phys. 13, 023019 (2011). 
[24| W.L. You, Y.W. Li, and S.J. Gu, Phys. Rev. E 76, 
022101 (2007). 

[25[ M. Cozzini, R. lonicioiu, P. Zanardi, Phys. Rev. B 76, 
104420 (2007). 

[26[ P. Buonsante, A. Vezzani, Phys. Rev. Lett. 98, 110601 

(2007) . 

[27[ G. Vidal, Phys. Rev. Lett. 101, 110501 (2008). 
[28| G. Evenbly and G. Vidal, Phys. Rev. B 79, 144108 
(2009). 

[29[ See supplementary material for technical details regard- 
ing implementation of MERA. 



Dipolar Molecules in Optical Lattices 
Supplementary Material 



Tomasz Sowinski, Omjyoti Dutta. Philipp Hauke. Luca Tagliacozzo, Maciej Lewenstein 



Calculation of hopping terms T and P 

Here we describe the procedure to calculate the terms 
in the modified Hubbard model in Eq. (1). First we 
find the lowest Bloch band for a single-particle mov- 
ing in the potential V\att{r) = Vq [sin^ + sin^ ^y] + 
From that, we construct the Wannier func- 
tions Wf{x,y)e~'^^ localized at site i By expanding 
the field operator in the Wannier basis, we derive the 
parameters for the Hubbard model. In particular, the 
integrals used to calculate the correlated hopping term T 
and the pair-hopping term P are: 

T = J y"d'Vd3r'yvf(x,?/)e-'=^'v(r-r') 

P=IJ J d^rd^r'W.ix, y)Wj [x, y)e-""' 
X V{r - r')W,ix', y')Wjix', yOe"""" 

(SI) 

Description of MERA 

MERA is a variational method that consists in pos- 
tulating a tensor-network structure for the low-energy 
states of Hamiltonian. The tensor network T is i) built 
from elementary tensors belonging to two different fami- 
lies, isometries Ii and disentanglers Vi that are isometric, 

laj = I: ■D.vj = I ; (S2) 

and ii) has a layered structure T ~ Yii Ti, such that 
each layer Ti performs an ER transformation (2I, |3| from 
a lattice Ci with lattice spacing bi to a lattice with 
spacing 6^+1 = nbi. Property ii) is at the origin of the 
ability of the MERA ansatz to describe infinite critical 
states with finite computational resources. This is the 
advantage of the MERA with respect to more traditional 
methods for studying ID chains such as, e.g., DMRG. 
Symmetries of the Hamiltonian can be encoded in the 
structure of the tensors. For example, in order to en- 
code translational invariant states of chains with periodic 
boundary conditions, we use inside each layer the same 
isometry and disentangler as many times as required to 
complete the ER transformation from the lattice Ci to 
the lattice Ci+i. When all the isometries and disenta- 
glers inside a given layer are chosen to be the same, the 
factor n not only characterize the blocking factor of the 
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FIG. SI. Tensors Ii (isometries) and Vi (disentanlgers) 
are represented by circles with trailing legs representing their 
indices. Lines connecting two tensors represent tensor con- 
tractions over the involved indices, i) The tensors are chosen 
such as to fuUfill the isometry constraints defined in Eq. (jS2|l . 
ii) A layer Ti of the 4 1 MERA tensor network T that 
maps operators and states defined on a lattice d with lattice 
spacing hi to operators and states defined on a lattice d+i 
with lattice spacing 46i Q. 



ER procedure (we talk about n to 1 MERA) but it also 
defines the size of the unit cell of the state. 

In the model we are considering the presence of CB 
patterns in some parts of the phase diagram extracted 
from ED suggests that we need an ansatz that can natu- 
rally encode at least a unit cell of two sites. This can be 
accomplished hy a 2 to 1 MERA, i.e., by blocking two 
sites into one at each step of the ER procedure. However, 
this MERA is computationally more expensive than the 3 
to 1 MERA. Unfortunately, the translationally invariant 
3 to 1 MERA does not easily accomodate a CB pattern, 
whence we choose a, 4 to 1 MERA that both naturally 
accomodates the two-site unit cell of a CB phase and re- 
duces the computational cost of the 2 to 1 MERA [Hj. 
In Fig. ISf f a) ii), we show a layer of the TN structure 
for the 4 to 1 MERA that we have used. MERA has a 
refinement parameter m larger values of which provide 
more accurate results but imply larger simulation time, 
since the complexity of the algorithm is ©(m^) in mem- 
ory and O(m^) in number of operations per iteration jH, 
and modest values of m such as m = 8 are often enough 
to get a correct qualitative picture of the model. 
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